### Alizade, Dancygier, Ditlmann
### "National Penalties Reversed"
### Replication Code 
### Figure A2
### For questions, contact jalizade@princeton.edu

# setup
rm(list = ls())
setwd("C:/Users/Jey/Dropbox/WZB/NaturalizationExperiment/Submission/JOP/replication_JOP/data")
library(foreign)
library(equivalence)
library(ggplot2)
dat <- read.dta("data_experimental.dta")
datmun <- read.dta("data_covars_mun.dta")



# function to run equivalence test
eqtest <- function(df, treat, cov) {
  tost(x=df[treat==0, cov], y=df[treat==1, cov], alpha = 0.05, var.equal = F, conf.level = 0.95, epsilon = 0.2*sd(df[,cov], na.rm =T))
}


# convert treatment variable from second experiment to factor
dat$e2_treat[dat$e2_treat=="NA"] <- NA
dat$e2_treat <- factor(dat$e2_treat, levels=c("Vote Intention (Canadian)", "Control (Turkish)", "Vote Intention (Turkish)", "Integration Problems (Turkish)"))

### Figure A2a ###

# binarize if necessary
datmun$med_pop <- ifelse(datmun$pop >= median(datmun$pop, na.rm = T), 1, 0)
datmun$med_dens <- ifelse(datmun$pop_density >= median(datmun$pop_density, na.rm = T), 1, 0)
datmun$e1_med_council <- ifelse(datmun$e1_n_council >= median(datmun$e1_n_council, na.rm = T), 1, 0)
dat$e1_age_dummy <- ifelse(dat$e1_agepol == 3, 1, 0)

# scale to 0-1 if necessary
datmun[,c("pct_forcit_pop", "pct_coll_qual", "pct_welfare", "e1_pct_left", "e1_pct_bigpart")] <-  datmun[,c("pct_forcit_pop", "pct_coll_qual", "pct_welfare", "e1_pct_left", "e1_pct_bigpart")] / 100

# subset first experiment
datmun1 <- datmun[datmun$e1_treated==1,]
dat1 <- dat[dat$e1_treated==1,]

# run tests
tests_a <- list(eqtest(df=datmun1, treat=datmun1$e1_treat_turkish, cov = "pct_welfare"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_turkish, cov = "e1_pct_left"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_turkish, cov = "med_pop"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_turkish, cov = "e1_pct_bigpart"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_turkish, cov = "e1_med_council"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_turkish, cov = "pct_coll_qual"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_turkish, cov = "med_dens"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_turkish, cov = "pct_forcit_pop"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_turkish, cov = "east"),
           eqtest(df=dat1, treat=dat1$e1_treat_turkish, cov = "e1_age_dummy"),
           eqtest(df=dat1, treat=dat1$e1_treat_turkish, cov = "e1_fempol"),
           eqtest(df=dat1, treat=dat1$e1_treat_turkish, cov = "e1_cfpol"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_dual, cov = "pct_welfare"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_dual, cov = "e1_pct_left"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_dual, cov = "med_pop"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_dual, cov = "e1_pct_bigpart"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_dual, cov = "e1_med_council"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_dual, cov = "pct_coll_qual"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_dual, cov = "med_dens"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_dual, cov = "pct_forcit_pop"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_dual, cov = "east"),
           eqtest(df=dat1, treat=dat1$e1_treat_dual, cov = "e1_age_dummy"),
           eqtest(df=dat1, treat=dat1$e1_treat_dual, cov = "e1_fempol"),
           eqtest(df=dat1, treat=dat1$e1_treat_dual, cov = "e1_cfpol"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_benefits, cov = "pct_welfare"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_benefits, cov = "e1_pct_left"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_benefits, cov = "med_pop"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_benefits, cov = "e1_pct_bigpart"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_benefits, cov = "e1_med_council"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_benefits, cov = "pct_coll_qual"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_benefits, cov = "med_dens"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_benefits, cov = "pct_forcit_pop"),
           eqtest(df=datmun1, treat=datmun1$e1_treat_benefits, cov = "east"),
           eqtest(df=dat1, treat=dat1$e1_treat_benefits, cov = "e1_age_dummy"),
           eqtest(df=dat1, treat=dat1$e1_treat_benefits, cov = "e1_fempol"),
           eqtest(df=dat1, treat=dat1$e1_treat_benefits, cov = "e1_cfpol"))

# differences in means
diffs_a <- sapply(tests_a, function(x) x$estimate[2]-x$estimate[1])

# results of tests
results_a <- sapply(tests_a, function(x) x$result)
results_a[results_a=="rejected"] <- "Rejected"

# create data frame with mean differences and results
df_a <- data.frame(diff=diffs_a, result=results_a)
df_a$treat <- rep(c("Turkish", "Dual", "Benefits"), each=12)
df_a$cov <- rep(c("Social Welfare Recipients", "Center-Left Parties", "Large Population", "Large Parties", 
                  "Large Council", "College Qualification", "High Population Density", "Foreign Population", 
                  "East German", "Politician 66+", "Female Politician", "Contact Form"), 3)
df_a$cov <- factor(df_a$cov, levels=rev(c("Social Welfare Recipients", "Politician 66+", "Center-Left Parties",
                                          "Large Population", "Large Parties", "Large Council", "College Qualification",
                                          "High Population Density", "Foreign Population", "Female Politician", 
                                          "East German", "Contact Form")))

# plot
fig_a <- ggplot(data=df_a, aes(x = cov, y= diff, group = treat, shape = treat, fill = result, color = result))
fig_a <- fig_a + geom_hline(yintercept=0, linetype="dashed")
fig_a <- fig_a + geom_point(position = position_dodge(width = 0.5))
fig_a <- fig_a + xlab("") + ylab("Difference in Means between Treatment Conditions")
fig_a <- fig_a + scale_y_continuous(limits = c(-0.05, 0.05), breaks = seq(-0.04, 0.04, 0.02))
fig_a <- fig_a + coord_flip()
fig_a <- fig_a + scale_color_manual(values ="black", labels = "Rejected", name = "Null Hypothesis\nof Difference", drop=FALSE)
fig_a <- fig_a + scale_fill_manual(values = "black", labels = "Rejected", name = "Null Hypothesis\nof Difference", drop=FALSE)
fig_a <- fig_a + scale_shape_manual(values = c(21, 22, 23), name = "Treatment")
fig_a <- fig_a + theme(legend.title=element_text(size=10, face = "bold"))
fig_a <- fig_a + guides(shape = guide_legend(reverse =T, order = 1))
fig_a
#ggsave("fgA2a.pdf", fig_a, width = 7, height=5)


### Figure A2b ###

# binarize if necessary
datmun$e2_med_council <- ifelse(datmun$e2_n_council >= median(datmun$e2_n_council, na.rm = T), 1, 0)
dat$e2_age_dummy <- ifelse(dat$e2_agepol == 3, 1, 0)

# scale to 0-1 if necessary
datmun[,c("e2_pct_left", "e2_pct_bigpart")] <-  datmun[,c("e2_pct_left", "e2_pct_bigpart")] / 100

# subset first experiment
datmun2 <- datmun[datmun$e2_treated==1,]
dat2 <- dat[dat$e2_treated==1,]

# construct treatment variables
datmun2$e2_canvote <- ifelse(datmun2$e2_treat=="Vote Intention (Canadian)", 1, ifelse(datmun2$e2_treat=="Vote Intention (Turkish)", 0, NA))
datmun2$e2_turkcontrol <- ifelse(datmun2$e2_treat=="Control (Turkish)", 1, ifelse(datmun2$e2_treat=="Vote Intention (Turkish)", 0, NA))
datmun2$e2_turkint <- ifelse(datmun2$e2_treat=="Integration Problems (Turkish)", 1, ifelse(datmun2$e2_treat=="Vote Intention (Turkish)", 0, NA))

dat2$e2_canvote <- ifelse(dat2$e2_treat=="Vote Intention (Canadian)", 1, ifelse(dat2$e2_treat=="Vote Intention (Turkish)", 0, NA))
dat2$e2_turkcontrol <- ifelse(dat2$e2_treat=="Control (Turkish)", 1, ifelse(dat2$e2_treat=="Vote Intention (Turkish)", 0, NA))
dat2$e2_turkint <- ifelse(dat2$e2_treat=="Integration Problems (Turkish)", 1, ifelse(dat2$e2_treat=="Vote Intention (Turkish)", 0, NA))

# run tests
tests_b <- list(eqtest(df=datmun2, treat=datmun2$e2_canvote, cov = "pct_welfare"),
                eqtest(df=datmun2, treat=datmun2$e2_canvote, cov = "e2_pct_left"),
                eqtest(df=datmun2, treat=datmun2$e2_canvote, cov = "med_pop"),
                eqtest(df=datmun2, treat=datmun2$e2_canvote, cov = "e2_pct_bigpart"),
                eqtest(df=datmun2, treat=datmun2$e2_canvote, cov = "e2_med_council"),
                eqtest(df=datmun2, treat=datmun2$e2_canvote, cov = "pct_coll_qual"),
                eqtest(df=datmun2, treat=datmun2$e2_canvote, cov = "med_dens"),
                eqtest(df=datmun2, treat=datmun2$e2_canvote, cov = "pct_forcit_pop"),
                eqtest(df=datmun2, treat=datmun2$e2_canvote, cov = "east"),
                eqtest(df=dat2, treat=dat2$e2_canvote, cov = "e2_age_dummy"),
                eqtest(df=dat2, treat=dat2$e2_canvote, cov = "e2_fempol"),
                eqtest(df=dat2, treat=dat2$e2_canvote, cov = "e2_cfpol"),
                eqtest(df=datmun2, treat=datmun2$e2_turkcontrol, cov = "pct_welfare"),
                eqtest(df=datmun2, treat=datmun2$e2_turkcontrol, cov = "e2_pct_left"),
                eqtest(df=datmun2, treat=datmun2$e2_turkcontrol, cov = "med_pop"),
                eqtest(df=datmun2, treat=datmun2$e2_turkcontrol, cov = "e2_pct_bigpart"),
                eqtest(df=datmun2, treat=datmun2$e2_turkcontrol, cov = "e2_med_council"),
                eqtest(df=datmun2, treat=datmun2$e2_turkcontrol, cov = "pct_coll_qual"),
                eqtest(df=datmun2, treat=datmun2$e2_turkcontrol, cov = "med_dens"),
                eqtest(df=datmun2, treat=datmun2$e2_turkcontrol, cov = "pct_forcit_pop"),
                eqtest(df=datmun2, treat=datmun2$e2_turkcontrol, cov = "east"),
                eqtest(df=dat2, treat=dat2$e2_turkcontrol, cov = "e2_age_dummy"),
                eqtest(df=dat2, treat=dat2$e2_turkcontrol, cov = "e2_fempol"),
                eqtest(df=dat2, treat=dat2$e2_turkcontrol, cov = "e2_cfpol"),
                eqtest(df=datmun2, treat=datmun2$e2_turkint, cov = "pct_welfare"),
                eqtest(df=datmun2, treat=datmun2$e2_turkint, cov = "e2_pct_left"),
                eqtest(df=datmun2, treat=datmun2$e2_turkint, cov = "med_pop"),
                eqtest(df=datmun2, treat=datmun2$e2_turkint, cov = "e2_pct_bigpart"),
                eqtest(df=datmun2, treat=datmun2$e2_turkint, cov = "e2_med_council"),
                eqtest(df=datmun2, treat=datmun2$e2_turkint, cov = "pct_coll_qual"),
                eqtest(df=datmun2, treat=datmun2$e2_turkint, cov = "med_dens"),
                eqtest(df=datmun2, treat=datmun2$e2_turkint, cov = "pct_forcit_pop"),
                eqtest(df=datmun2, treat=datmun2$e2_turkint, cov = "east"),
                eqtest(df=dat2, treat=dat2$e2_turkint, cov = "e2_age_dummy"),
                eqtest(df=dat2, treat=dat2$e2_turkint, cov = "e2_fempol"),
                eqtest(df=dat2, treat=dat2$e2_turkint, cov = "e2_cfpol"))


# differences in means
diffs_b <- sapply(tests_b, function(x) x$estimate[2]-x$estimate[1])

# results of tests
results_b <- sapply(tests_b, function(x) x$result)
results_b[results_b=="rejected"] <- "Rejected"

# create data frame with mean differences and results
df_b <- data.frame(diff=diffs_b, result=results_b)
df_b$treat <- rep(c("Vote Intention (Canadian)", "Control (Turkish)", "Integration\nProblems (Turkish)"), each=12)
df_b$treat <- factor(df_b$treat, levels = rev(c("Control (Turkish)", "Vote Intention (Canadian)", "Integration\nProblems (Turkish)")))
df_b$cov <- rep(c("Social Welfare Recipients", "Center-Left Parties", "Large Population", "Large Parties", 
                  "Large Council", "College Qualification", "High Population Density", "Foreign Population", 
                  "East German", "Politician 66+", "Female Politician", "Contact Form"), 3)
df_b$cov <- factor(df_b$cov, levels=rev(c("Social Welfare Recipients", "Politician 66+", "Center-Left Parties",
                                          "Large Population", "Large Parties", "Large Council", "College Qualification",
                                          "High Population Density", "Foreign Population", "Female Politician", 
                                          "East German", "Contact Form")))


# plot
fig_b <- ggplot(data=df_b, aes(x = cov, y= diff, group = treat, shape = treat, fill = result, color = result))
fig_b <- fig_b + geom_hline(yintercept=0, linetype="dashed")
fig_b <- fig_b + geom_point(position = position_dodge(width = 0.5))
fig_b <- fig_b + xlab("") + ylab("Difference in Means between Treatment Conditions")
fig_b <- fig_b + scale_y_continuous(limits = c(-0.05, 0.05), breaks = seq(-0.04, 0.04, 0.02))
fig_b <- fig_b + coord_flip()
fig_b <- fig_b + scale_color_manual(values =c("red", "black"), labels = c("Not Rejected", "Rejected"), name = "Null Hypothesis\nof Difference", drop=FALSE)
fig_b <- fig_b + scale_fill_manual(values =c("red", "black"), labels = c("Not Rejected", "Rejected"), name = "Null Hypothesis\nof Difference", drop=FALSE)
fig_b <- fig_b + scale_shape_manual(values = c(21, 22, 23), name = "Treatment")
fig_b <- fig_b + theme(legend.title=element_text(size=10, face = "bold"))
fig_b <- fig_b + guides(shape = guide_legend(reverse =T, order = 1))
fig_b
#ggsave("fgA2b.pdf", fig_b, width = 7, height=5)
